### =========================================================================
### Preparations
### =========================================================================

### Indicate the working path
wkpath <- ""
setwd(wkpath)


### Please indicate the path to the environmental layers
SDM_env_path <- ""
### The 10 environmental layers can be downloaded from:
### https://www.envidat.ch/dataset/sdm-env-layers-gdplants


### Please set the folder with the occurrence data. 
### Data should contain two columns "x" and "y" to give the longitude and latitude of occurrences
input_dir <- "2.2_cleaning_cc/cleaning_cc_occurrences"

### Please set the folder for the output maps
output_dir <- "3.2_mapping_SDM/0_parameter_optimization"

### Please load polygon mapping functions from where you save these scripts
cmon.files=list.files("GDPlants/functions/sdm_functions",full.names = T)
sapply(cmon.files,source)

### Load Kew data
### Please contact Kew (https://powo.science.kew.org/) for distribution data and map
### Or you can use your own reference database
kewdata <- read.table("./Kew_data_world_plant/Plant_Data_Kew.txt", header=T)
kewmap <- readOGR(dsn = "./Kew_data_world_plant/KewMap/level3.shp")


### =========================================================================
### Preparations of packages, functions and data ####
### =========================================================================

# install packages globally

# set a local mirror server
options(repos=structure(c(CRAN="https://stat.ethz.ch/CRAN/")))

# Packages to load
packages <- c("raster","rgdal","maptools", "tools","dismo","cluster","class",
              "gam","gbm","randomForest","ROCR","parallel","spatstat","rgeos")

# Install missing ones and load all packages
for (p in packages) {
  if(!library(package = p, logical.return = TRUE, character.only = TRUE)){
    install.packages(p)
    library(package = p, character.only = TRUE)
  } else {   
    library(package = p, character.only = TRUE) 
  }
}


### Prepare Kew data

kewdata_native <- kewdata[kewdata$Introduced == 0,]
kewdata_family_table <- na.omit(kewdata_native)

kewdata_family_table$name <- paste(kewdata_family_table$Genus,kewdata_family_table$Species)

# =========================================================================


### Extract species families and names ####
allfiles <- list.files(input_dir, full.names = F, pattern = ".\\csv")

allfilenames <- sapply(allfiles, function(x)   {strsplit(x, ".\\csv")[[1]][1]})

spfamily <- sapply(allfilenames, function(x)   {strsplit(x, "_")[[1]][1]})
Spname <- sapply(allfilenames, function(x)   {strsplit(x, "_")[[1]][2]})


### =========================================================================
### Prepare validation dataset based on real distribution from Kew ####
### NOTICE: if you have generated the validation points before, you can skip this step
### =========================================================================

if (!dir.exists(file.path(output_dir, "validation_results"))){
  dir.create(file.path(output_dir, "validation_results"), recursive = T)
}
if (!dir.exists(file.path(output_dir, "validation_points_from_kew"))){
  dir.create(file.path(output_dir, "validation_points_from_kew"), recursive = T)
}


for (id in 1:length(Spname)){
  
  spi_name <- Spname[id]
  
  if (spi_name %in% kewdata_family_table$name){
  
  kewsprange <- kewdata_family_table[kewdata_family_table$name == spi_name, ]
  kewsprangemap <- kewmap[which(kewmap@data$LEVEL3_COD %in% unique(kewsprange$Area_code_L3)),]
  
  kewsprangemap_abs <- kewmap[which(!kewmap@data$LEVEL3_COD %in% unique(kewsprange$Area_code_L3)),]
  
  
  # generate presences and absences based on real distribution
  validation_prs <- spsample(kewsprangemap,2000, "random")
  validation_abs <- spsample(kewsprangemap_abs,2000, "random")
  validation_points <- as.data.frame(rbind(validation_prs,validation_abs))
  validation_points$Presence <- rep(c(1,0), each=2000)
  
  write.csv(validation_points, file.path(output_dir, "validation_points_from_kew", paste0(spfamily[id],"_",spi_name,".csv")), row.names=FALSE)
  
  } else {
    print(paste(spi_name, "is not recorded by Kew."))
  }
}


### =========================================================================
### Load environmental variables ####
### =========================================================================

# Load environmental layers

env.stk= stack(file.path(SDM_env_path, 'CHELSA_bio10_01_preped.tif'),
               file.path(SDM_env_path, 'aridity_tranformed.tif'),
               file.path(SDM_env_path, 'ORCDRC_M_sl1_1km_ll_tranformed.tif'),
               file.path(SDM_env_path, 'CHELSA_FCF_1979-2013_V1.2_preped.tif'),
               file.path(SDM_env_path, 'CHELSA_bio10_17_tranformed.tif'),
               file.path(SDM_env_path, 'PHIHOX_M_sl1_1km_ll_preped.tif'),
               file.path(SDM_env_path, 'CHELSA_bio10_02_preped.tif'),
               file.path(SDM_env_path, 'CHELSA_bio10_15_preped.tif'),
               file.path(SDM_env_path, 'CLYPPT_M_sl1_1km_ll_preped.tif'))

pred_sdm <- c('Mean_Temp', 'Aridity', 'Organic_content_soil',"Frost_Change_Frequency","Prec_Driest_Quarter",
              "Soil_pH","Mean_Diurnal_Range","Prec_Seasonality","Clay_content_soil") 
names(env.stk)=pred_sdm

proj <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs +towgs84=0,0,0")


### =========================================================================
### Loop over species and create maps ####
### =========================================================================

PAmethod <- c('target.group','target.group.allsp','geographic', 'density', 'random',  'geo.strat', 'env.strat', 'env.semi.strat') #


#### Start testing ####

for (id in 1:length(Spname)){
  
  # Get name
  env.stk_i=env.stk
  spi_name <- Spname[id]
  print(paste(id,spi_name))
  
  #load the presence
  spp_occ <- read.csv(list.files(input_dir, pattern = paste0(spi_name, ".csv"), full.names = T))
  
  if (nrow(spp_occ) >= 20){
    
    spp_occ <- spp_occ[,c("x","y")]
    spp_occ <- SpatialPoints(spp_occ, proj4string=proj)
    
    for (PAmn in 1:length(PAmethod)){
      
      PAm = PAmethod[PAmn]
      
      if (PAm == 'target.group.allsp'){
        
        pseu.abs_i <- wsl.samplePseuAbs(type="target.group",
                                        n=10000, #min 1000
                                        env.stack=env.stk_i,
                                        add.strat=0.2,
                                        env.strat_path = file.path(output_dir, "envpath"),
                                        pres = spp_occ,
                                        template_dir= file.path(output_dir, "template_dir"),
                                        target.group_dir = file.path(input_dir, unique(spfamily)),
                                        taxon=spi_name,
                                        force_spat_thin = "both",
                                        geores_fact=10,
                                        limdist=NA)
        print(paste("Create pseudoabsences for", spi_name, "with target.group strategy (with all reference species)."))
        
      } else { # For the other PA strategies
        
        pseu.abs_i <- wsl.samplePseuAbs(type=PAm,
                                        n=10000, #min 1000
                                        env.stack=env.stk_i,
                                        add.strat=0.2,
                                        env.strat_path = file.path(output_dir, "envpath"),
                                        pres = spp_occ,
                                        template_dir = file.path(output_dir, "template_dir"),
                                        target.group_dir = file.path(input_dir,spfamily[id]),
                                        taxon=spi_name,
                                        force_spat_thin = "both",
                                        geores_fact=10,
                                        limdist=NA)
      }
      
      print(paste("Create pseudoabsences for", spi_name, "with", PAm, "strategy."))
      
      ### =========================================================================
      ### Decide on predictor number depending on available presence observations ####
      ### =========================================================================

      ssize=floor(length(which(pseu.abs_i@pa==1))/10)
      
      # If there are not at least 20 presence observations, do nothing
      if(ssize<2){
        next
        # If there are less than 90 observations reduce predictor set so that
        # at least 10 observations are available per predictor
      } else if(ssize<9){
        env.stk_i=env.stk[[1:ssize]]
        pseu.abs_i@env_vars=pseu.abs_i@env_vars[, 1:ssize,drop=F]
      }
      
      print("new pseu.abs_i created")
      
      ### =========================================================================
      ### Prepare data for modelling ####
      ### =========================================================================

      wip=which(pseu.abs_i@pa==1)
      
      blkp=sample(rep(1:3,each=ceiling(length(wip)/3)),ceiling(length(wip)/3)*3)[1:length(wip)]
      wia=which(pseu.abs_i@pa==0)
      blka=sample(1:3,size = length(wia),replace = T)
      
      blks=rep(NA,length(pseu.abs_i@pa))
      blks[wip]=blkp
      blks[wia]=blka
      
      ### =========================================================================
      ### Define SDMs ####
      ### =========================================================================
      
      ### Define formulas
      form.glm.s=as.formula(paste("Presence~",paste(names(env.stk_i),collapse="+"))) # simple
      form.glm.s2=as.formula(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",2)"),collapse="+"))) # intermediate
      
      if(nlayers(env.stk_i)>1){
        cmbs<-combn(names(env.stk_i),2)
        pst<-apply(cmbs,2,paste,collapse=":")
        int.part=paste(pst,collapse="+")
        form.glm.i=as.formula(paste(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",3)"),collapse="+")),int.part,sep="+")) # complex
      } else {
        form.glm.i=as.formula(paste(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",3)"),collapse="+")))) # complex
      }
      
      if(nlayers(env.stk_i)>1){
        cmbs<-combn(names(env.stk_i),2)
        pst<-apply(cmbs,2,paste,collapse=":")
        int.part=paste(pst,collapse="+")
        form.glm.i2=as.formula(paste(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",4)"),collapse="+")),int.part,sep="+")) # complex
      } else {
        form.glm.i2=as.formula(paste(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",4)"),collapse="+")))) # complex
      }
      
      if(nlayers(env.stk_i)>1){
        cmbs<-combn(names(env.stk_i),2)
        pst<-apply(cmbs,2,paste,collapse=":")
        int.part=paste(pst,collapse="+")
        form.glm.c=as.formula(paste(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",5)"),collapse="+")),int.part,sep="+")) # complex
      } else {
        form.glm.c=as.formula(paste(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",5)"),collapse="+")))) # complex
      }
      
      if(nlayers(env.stk_i)>1){
        cmbs<-combn(names(env.stk_i),2)
        pst<-apply(cmbs,2,paste,collapse=":")
        int.part=paste(pst,collapse="+")
        form.glm.c2=as.formula(paste(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",6)"),collapse="+")),int.part,sep="+")) # complex
      } else {
        form.glm.c2=as.formula(paste(paste("Presence~",paste(paste0("poly(",names(env.stk_i),",6)"),collapse="+")))) # complex
      }
      
      form.gam.s=as.formula(paste("Presence~",paste(paste0("s(",names(env.stk_i),",df=1)"),collapse="+"))) # simple
      form.gam.s2=as.formula(paste("Presence~",paste(paste0("s(",names(env.stk_i),",df=2)"),collapse="+"))) # simple
      form.gam.i=as.formula(paste("Presence~",paste(paste0("s(",names(env.stk_i),",df=3)"),collapse="+"))) # intermediate
      form.gam.i2=as.formula(paste("Presence~",paste(paste0("s(",names(env.stk_i),",df=5)"),collapse="+"))) # simple
      form.gam.c=as.formula(paste("Presence~",paste(paste0("s(",names(env.stk_i),",df=10)"),collapse="+"))) # complex
      form.gam.c2=as.formula(paste("Presence~",paste(paste0("s(",names(env.stk_i),",df=15)"),collapse="+"))) # complex
      
      form.trees=Presence ~ .
      
      modinp=list(multi("glm",list(formula=form.glm.s,family="binomial"),"glm-simple",step=FALSE,weight=FALSE),
                  multi("glm",list(formula=form.glm.s2,family="binomial"),"glm-simple2",step=FALSE,weight=FALSE),
                  multi("glm",list(formula=form.glm.i,family="binomial"),"glm-interm",step=FALSE,weight=FALSE),
                  multi("glm",list(formula=form.glm.i2,family="binomial"),"glm-interm2",step=FALSE,weight=FALSE),
                  multi("glm",list(formula=form.glm.c,family="binomial"),"glm-complex",step=FALSE,weight=FALSE),
                  multi("glm",list(formula=form.glm.c2,family="binomial"),"glm-complex2",step=FALSE,weight=FALSE),
                  multi("gam",list(formula=form.gam.s,family="binomial"),"gam-simple",step=FALSE,weight=FALSE),
                  multi("gam",list(formula=form.gam.s2,family="binomial"),"gam-simple2",step=FALSE,weight=FALSE),
                  multi("gam",list(formula=form.gam.i,family="binomial"),"gam-interm",step=FALSE,weight=FALSE),
                  multi("gam",list(formula=form.gam.i2,family="binomial"),"gam-interm2",step=FALSE,weight=FALSE),
                  multi("gam",list(formula=form.gam.c,family="binomial"),"gam-complex",step=FALSE,weight=FALSE),
                  multi("gam",list(formula=form.gam.c2,family="binomial"),"gam-complex2",step=FALSE,weight=FALSE),
                  multi("gbm",list(formula=form.trees,
                                   distribution = "bernoulli",
                                   interaction.depth = 5,
                                   shrinkage=.005,
                                   n.trees = 100),"gbm-simple",weight=FALSE),
                  multi("gbm",list(formula=form.trees,
                                   distribution = "bernoulli",
                                   interaction.depth = 5,
                                   shrinkage=.005,
                                   n.trees = 200),"gbm-simple2",weight=FALSE),
                  multi("gbm",list(formula=form.trees,
                                   distribution = "bernoulli",
                                   interaction.depth = 5,
                                   shrinkage=.005,
                                   n.trees = 300),"gbm-interm",weight=FALSE),
                  multi("gbm",list(formula=form.trees,
                                   distribution = "bernoulli",
                                   interaction.depth = 5,
                                   shrinkage=.005,
                                   n.trees = 500),"gbm-interm2",weight=FALSE),
                  multi("gbm",list(formula=form.trees,
                                   distribution = "bernoulli",
                                   interaction.depth = 5,
                                   shrinkage=.005,
                                   n.trees = 1000),"gbm-complex",weight=FALSE),
                  multi("gbm",list(formula=form.trees,
                                   distribution = "bernoulli",
                                   interaction.depth = 5,
                                   shrinkage=.005,
                                   n.trees = 10000),"gbm-complex2",weight=FALSE),
                  multi("randomForest",list(formula=form.trees,ntree=500,nodesize=100),"rf-simple"),
                  multi("randomForest",list(formula=form.trees,ntree=500,nodesize=50),"rf-simple2"),
                  multi("randomForest",list(formula=form.trees,ntree=500,nodesize=30),"rf-interm"),
                  multi("randomForest",list(formula=form.trees,ntree=500,nodesize=10),"rf-interm2"),
                  multi("randomForest",list(formula=form.trees,ntree=500,nodesize=3),"rf-complex"),
                  multi("randomForest",list(formula=form.trees,ntree=500,nodesize=1),"rf-complex2"))
      
      ### =========================================================================
      ### fit models for prediction ####
      ### =========================================================================
      print(paste("Start model fitting for", spi_name))
      prmod=wsl.flex(x=pseu.abs_i,
                     replicatetype="none",
                     reps=1,
                     project="tree_map",
                     mod_args=modinp)
      
      ### =========================================================================
      ### Validate independently ####
      ### =========================================================================

      print(paste("Start validation for", spi_name, "with", PAm, "strategy."))
      
      validation_points <- read.csv(list.files(file.path(output_dir, spfamily[id]), pattern = spi_name), head=T)
      
      validation_env <- extract(env.stk_i, validation_points[,c("x","y")])
      tester <- na.omit(cbind(validation_points, validation_env))
      
      evalsValidate<-wsl.evaluate(prmod,crit="maxTSS",prevalence_correction = T, tester = tester)
      
      smev_val<-summary(evalsValidate)
      
      write.csv(smev_val, 
                file.path(output_dir, "validation_results", paste0(PAm,"_", spfamily[id], "_" , spi_name, ".csv")),
                row.names = TRUE)
      
      cat(paste("Validation result for", spi_name, "with", PAm, "strategy is saved in:\n",
                file.path(output_dir, "validation_results",paste0(PAm,"_", spfamily[id], "_" , spi_name, ".csv"))
      ))

    }
    print(paste0("COMPLETED:", id, "/", tohere))#round(id/length(tohere),2)*100, "% ]"))
    tmp <- tempfile()
    do.call(file.remove, list(list.files("tmp", full.names = TRUE)))
    
    
  }
}